SPACE VARYING GAMUT MAPPING 



1 Technical Field 

2 The technical field is color image processing using gamut correction. 

3 Background 

4 Gamut mapping is a method used to modify a representation of a color image to 

5 fit into a constrained color space of a given rendering medium. A laser-jet color printer 

6 that attempts to reproduce a color image on regular paper would have to map the 

7 photographed picture colors in a given color range, also known as the image "color 

8 gamut," into the given printer/page color gamut. Conventional gamut mapping methods 

9 involve a pixel by pixel mapping (usually a pre-defined look-up table) and ignore the 

10 spatial color configuration. More recently, spatial dependent approaches were proposed 

1 1 for gamut mapping. However, these solutions are either based on heurestic assumptions 

12 or involve a high computational cost. An example of such a method is discussed in 

13 "Color Gamut Mapping Based on a Perceptual Image Difference Measure," S. Nakauchi, 

14 et al., Color Research and Application . Vol. 24, pp. 280-291 (1999). 

1 5 Another method relies on preserving the magnitude of gradients in an original 

1 6 image, while projecting the gradients onto the target gamut as a constraint. This multi- 

1 7 scale property is achieved by sampling the image around each pixel with exponentially 

1 8 increasing sampling intervals while the sampling is done along vertical and horizontal 

1 9 axes. The method involves modulating an L 2 measure for image difference by human 

20 contrast sensitivity functions. The method uses a model in which the contrast sensitivity 

2 1 function is a linear combination of three spatial band-pass filters H } , H 2 , H 3 , given in the 

22 spatial-frequency domain {oxh 1 ,h 2 ,h 3 , as their corresponding spatial filters), as shown in 

23 Figure 1. 

24 For gamut mapping of the image u 0 in the CIELAB (which stands for Commission 

25 International de I'Eclairage (CIE) with L, A, and B (more commonly L, a, b) respectively 

26 representing converted signals for cyan (C), magenta (M) and yellow (Y)) space, the 

27 method minimizes the functional 
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1 subject to {u L , u, u b ) e 3. 

2 In Equation (1), is the filter corresponding to the spectral channel c e {£, a, &} the i 

3 e {1, 2, 3} 'contrast sensitivity' mod, Q, is the image domain, and <9 is the target gamut. 

4 Note that a total of nine filters are involved, three for each spectral channel and a 

5 total of three spectral channels. 

6 The filters H c t are modeled by shifted Gaussians. is not shifted, and thus, ti; 

7 is also Gaussian, while and h% are Gaussians modulated by two sine functions with 

8 different periods. A graphical analysis of and h% as shown in Figure 2 reveals that 

9 and h% approximate the derivative operator at different scales. These two gradient 

10 approximation operators are denoted by and V^ 2 . Note that any band pass filter can 

11 be considered as a version of a derivative operator. Furthermore, one possible extension 

12 of the ID derivative to 2D is the gradient. Thus, the minimization of the Equation (1) 

13 functional is similar to minimizing the following functional for each channel separately. 

^ *(u-u 0 ) 2 +|v^(«-tt 0 } 2 +|V< 2 ( M - Mo ) 2 rfQ (2) 

a 

14 Roughly speaking, the first term corresponds to the S-CIELAB perceptual 

1 5 measure, while the next two terms capture the need for matching the image variations at 

1 6 two selected scales that were determined by human perception models. One technical 

17 difficulty of the spatial filters corresponding to Equation (1) is their large numerical 

1 8 support. 

1 9 Another simple spatial-spectral measure for human color perception was proposed 

20 in "A Spatial Extension of CIELAB for Digital Color Image Reproduction," Zhang et al., 

21 Proceedings of the SID Symposium , Vol. 27, pp. 731-34, (1996). The 'S-CIELAB' 

22 defines a spatial -spectral measure for human color perception by a composition of spatial 

23 band-pass linear filters in the opponent color space followed by the CIELAB Euclidean 

24 perceptual color. 
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1 Summary 

2 Gamut mapping is a method to modify representation of a color image to fit into a 

3 constrained color space of a given rendering medium. A laser-jet color printer that 

4 attempts to reproduce a color image on regular paper would have to map the 

5 photographed picture colors in a given color range, also known as the image "color 

6 gamut," into the given printer/page color gamut. 

7 A method an apparatus uses a variational approach for space dependent gamut 

8 mapping. The method presents a new measure for the problem. A solution to the 

9 formulation according to the method is unique if the gamut of the target device is known. 

10 A quadratic programming efficient numerical solution provides real-time results. 

11 In an embodiment, the method comprises receiving an input image, converting 

12 color representations of an image pixel set to produce a corresponding electrical values 

13 set, applying a space varying algorithm to the electrical values set to produce a color- 

14 mapped value set, and reconverting the color-mapped value set to an output image. The 

1 5 space varying algorithm minimizes a variational problem represented by 

16 E(u) = J n (d 2 + a\ VD\ 2 )iQx , subject to u e & , where Q is a support of the image, 

17 a is a non-negative real number, D=g*(u-u 0 ), and g is a normalized Gaussian kernel with 

1 8 zero mean and a small variance a. 

1 9 The method may be implemented in a variety of image capture reproduction 

20 devices, including cameras and printers. The method may be embodied on a computer- 

21 readable medium, such as a CD-ROM, for example. 

22 Description of the Drawings 

23 The detailed description will refer to the following drawings, in which like 

24 numerals refer to like objects, and in which: 

25 Figure 1 is a qualitative description of filters modeling the human contrast 

26 sensitivity functions in the spectral domain; 

27 Figure 2 illustrates a shifted Gaussian; 

28 Figure 3 is a block diagram of an apparatus that uses a variational approach for 

29 gamut mapping; 

30 Figures 4-7 illustrate algorithms and processes used with the apparatus of Figure 

31 3; 

32 Figure 8 illustrates an example of a behavior of an algorithm used by the 
3 3 apparatus of Figure 3 ; 
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1 Figure 9 is a flowchart of an algorithm illustrating an alternative operation of the 

2 apparatus of Figure 3; and 

3 Figure 1 0 is a block diagram of a computer system that implements the algorithm 

4 of Figure 4. 

5 Detailed Description 

6 Gamut mapping is a method used to modify a representation of a color image to 

7 fit into a constrained color space of a given rendering medium. A laser-jet color printer 

8 that attempts to reproduce a color image on regular paper would have to map the 

9 photographed picture colors in a given color range, also known as the image "color 

10 gamut," into the given printer/page color gamut. Conventional gamut mapping methods 

1 1 involve a pixel by pixel mapping (usually a pre-defined look-up table) and ignore the 

12 spatial color configuration. More recently, spatial dependent approaches were proposed 

13 for gamut mapping. However, these solutions are either based on heurestic assumptions 

14 or involve a high computational cost. 

15 The gamut mapping problem is related to the Retinex problem of illumination 



16 compensation and dynamic range compression. The basic Retinex problem is: How to 

1 7 estimate the reflectance image from the given acquired image? A reasonable optical 

1 8 model of the acquired image S asserts that it is a multiplication of the reflectance R and 

19 the illumination L images. Where the reflectance image is a hypothetic image that would 

20 have been measured if every visible surface had been illuminated by a unit valued white 

21 illumination source, and the illumination image is the actual illumination shaded on 

22 surfaces in the scene. In the log domain: 



23 s = r + / 

24 where s, r, and / are the respective logarithms of S, R, and L. Since surface patches can 

25 not reflect more light than has been shaded on them, R < 1 => r < 0. Thus, in an ideal 

26 situation, image r < 0, which is perceptually similar to s. For the Retinex, an additional 

27 physically motivated constraint is provided, namely, that the illumination image I = s — r 

28 is smooth, i.e. the gradient | V/ 1 = | V(r -s) \ is small. But this is just another way to say 

29 that the features of r are similar to those of s, since the illumination is assumed not to 

30 create perceptual features in s. In the gamut mapping problem an image uo exists, and the 

3 1 problem is to estimate an image u s& , which is not only perceptually similar to u 0 , but 

32 also has similar perceptual features as uq. 
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1 A good measure of image deviation captures the perceptual difference between 

2 the initial, uq and the final, u, images. This is modeled by 

D = g*{u-u 0 ). (3) 

3 where g may be a normalized Gaussian kernel with zero mean and a small variance cr . 

4 This model is good for small deviations. However, for large deviations it should be 

5 elaborated to account for possible perceptual feature differences, which may be modeled 

6 by the difference of gradients, which due to linearity, turns out to be the gradient of 

7 Equation (3) 

VD = V[g*{u-u 0 )] = g*(Vu-Vu 0 ) (4) 

8 The proposed measure yields the functional 

E(u)= ^D 2 +a\VD\ 2 )dQ, (5) 

n 

9 subject to u<= 3, 

10 which should be minimized. 

1 1 Taking a first variation of Equation (5) with respect to u yields the Euler-Lagrange 

12 equation: 

= g * (adivfrg * <« - « 0 ))-**("- «o ))= 0> (6) 

ou 

13 Reformulating Equation (6) as a gradient descent flow for u, provides the following 

14 minimization scheme: 

s.t. ue3 

15 The functional (Equation (5)) and the resulting minimization scheme are both 

1 6 Euclidean invariant in the image plane. They are thus both translation and rotation 

17 invariant. As the parameter a goes to zero, the S-CIELAB model is approximated, while 

1 8 for effective a the result is a proper extension to the perceptual measures, with an 

19 efficient numerical implementation. 
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1 To provide a numerical solution, recall that an artificial time parameter t was 

2 added to the image u(x, y), which now reads u{x, y, i). A first step is to discretize the 

3 Eueler- Lagrange gradient descent equation, by first taking a simple forward explicit 

4 approximation for the / derivative, 

— = ag* AD -g* D 

r 

5 where r = dt and u n {x,y)& u(x,y;nr). 

6 Next, the space derivatives are solved. Let u n l} ~ u{ih,jh, nf), where uniform 

7 spatial spacings are assumed in the x and y directions of size h. Using central derivatives 

8 in space, 

U rr « Dm = -i±l ^ 



9 and the same for the y direction. Using the relation g * (g * u) = g x * g x * u, the 

10 algorithm next computes the kernels D2 = g x * g x + gy * gy = Dxg * Dxg + D y g * D y g. The 

1 1 explicit approximation reads: 

£>" = g*g*(u"-U Q ) 

L" = D 2 *^"-w 0 ) 
u y +1 = u "j + T ~ D i" ) 

12 subject to the constraint u" e 3. 

13 To speed up convergence, a standard coarse to fine pyramidal approach may be 

14 preferred. For example, an original image S may be composed of a set of 5 12 x 5 12 

15 pixels. The image S may then be averaged down (sub-sampled or decimated) by taking, 

16 for example, every 4 adjacent pixels and replacing the block of four pixels with a single 

1 7 pixel having a value that represents an average of the four pixels. This results in an image 

1 8 set that is 256 x 256 pixels. The sub-sampling can be continued, yielding a subsequent 

19 reduced pixel set, or layer, comprising 128 x 128 pixels. This process may continue, with 

20 each resulting layer being an increasingly coarser version of the original image S. The 

21 final layer will be a set of only 2x2 pixels, and likely will be blurred. The process of 
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1 sub-sampling or decimating to provide increasingly coarser details results in an image 

2 pyramid. The original image S represent the base of the image pyramid at 5 12 x 5 12 

3 pixels, followed by layers of 256 x 256 pixels, 128 x 128 pixels, 64 x 64 pixels, and so 

4 on. Using this approach of an image pyramid allows for rapid convergence of the 

5 solution to the functional of Equation (5). The process starts by using the coarsest layer 

6 of the image pyramid, and running a number of iterations until a solution is achieved. 

7 The solution is then used to initialize the next finer layer of the image pyramid, resulting 

8 in less computation in order to arrive at a solution. The process continues until the finest 

9 layer (e.g., the original, full resolution, image set for S) is reached. As an alternative to 

1 0 the image pyramid, a full multi-grid method may be used. 

1 1 The functional of Equation (5) has a Quadratic Programming (QP) form, since the 

12 penalty term is quadratic and the constraint is linear. If the set $ is convex, the overall 

13 problem is convex if and only if the Hessian of the functional Equation (5) is positive 

14 definite. In such a case, there is a unique local minimum that is also the global solution to 

1 5 the problem. In the above embodiment, the Hessian is given by g * (1 - a A) * g, which is 

16 positive definite for all a > 0. Thus, for a convex target gamut $ , there exists a unique 

1 7 solution. 

1 8 The above-described functional (Equation (5)) may be used in imaging devices to 

1 9 map a gamut of an image to the gamut of a device that attempts to reproduce the image. 

20 Figure 3 is a block diagram of an apparatus 100 that may be used for gamut mapping 

21 according to the functional of Equation (5). An image capture device 101 receives an 

22 input original image and produces an output image S 102, which is an electrical signal 

23 representing a colorimetric value image. For example, the capture device 101 may 

24 convert at least three colorimetric values of each pixel of the original image into 

25 corresponding electrical signals. The electrical signals may indicate the L, a, b values, for 

26 example. Other colorimetric values may be the XYZ tristimilus value and the L, U, V or 

27 device dependent RGB values. A gamut mapper 105 produces an in-gamut image S' 106. 

28 Finally, a rendering module 107 provides a rendered image S" 108. The rendering 

29 module 1 07 may be implemented as a color laser printer. The thus-generated image 5"' 

30 108 may represent a best-fit image, given gamut limitations of the device in which the 

3 1 image reproduction is to occur. 

32 The image S 1 02 may represent the image as sensed by the capture module 1 03 . 

33 The gamut mapper 1 05 applies an algorithm to extract and map the values of the image S 

34 102 into the gamut of the image reproduction device 107. In particular, the gamut 



HP 10005732 



-6- 



1 mapper 1 05 may apply an algorithm that solves the problem represented by Equation (5), 

2 thereby solving the gamut mapping problem and optimizing the output image S" 1 08. 

3 Figure 4 is a block diagram showing routines of an algorithm 120 that may be 

4 used to complete the gamut mapping function. In an embodiment, the algorithm 120 

5 provides a numerical solution to the functional of Equation (5). 

6 The algorithm 120 begins with input block 121, where an input to the algorithm 

7 120 is an image S of size [N, M], and the two parameters a and p\ 

8 In initialization block 123, a Gaussian pyramid of the image s is computed. The 



9 thus-constructed pyramid contains p resolution layers (from fine (1) to coarse (p)) with 

10 the current resolution layer (k) set to p. In block 125, T^ iterations of a gradient descent 

1 1 algorithm are applied to the resolution layer, until all resolution layers are checked, 

12 block 127. In block 129, the next resolution layer is updated. When all resolution layers 

1 3 are processed, the result is the final output of the algorithm 120. 

14 Figure 5 is a block diagram showing the initialization routine 123 in more detail. 

1 5 In block 1 3 1 , a counter is initialized with i = 1 . In block 132, the Gaussian pyramid is 

1 6 constructed by smoothing the image with a kernel, such as the kernel k PY R. 
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17 In block 133, the image is decimated by, for example, a 2:1 ratio. The process is repeated 

18 (block 134) and the counter is incremented (block 135) p times where p < lg 2 (min (M, 

19 N)). This process produces a sequence of images {S k } p k=l conventionally named the 

20 "Gaussian Pyramid of 5"." The image Si is the original image S, and S p is an image with 

21 the coarsest resolution for the Gaussian pyramid. 

22 In block 137, the process sets k=p, i.e., starts at the coarsest resolution layer k, 

23 and sets the initial condition L 0 = max {S p } . 

24 Figure 6 is a block diagram of an embodiment of the processing main routine 125. 

25 The routine 125 starts at block 141 , where the new resolution layer is initialized and Tk 

26 and a.k are set (e.g., T k = K*To). 
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Then, for 7 = 1, .., 7k. 

In block 143, the routine 125 calculates the gradient: 
G = A(u -u 0 )+a k (u-u 0 ) 



2 where Ax is the convolution of each of the color planes of x with K L ap: 
and a k is, for example a k = a *2 2(k_1) 



0 1 0 

1 -2 1 
0 1 0 



In block 145, the routine 125 calculates jumsd- 



--ZG 2 / 



In block 147, the routine 125 completes the gradient descent iteration: 

Where u 0 is a constant; for example u o =0.8. 
In block 148, the result is projected onto the constraints: 

L,=Proj,(zJ 
Where Proj 9 (x) is a projection of x into the gamut 3. 

In block 149, if j * T K , processing returns to block 143. Otherwise, the routine 
125 ends. 

Returning to Figure 4, in block 127, if k > 1, the result L n is up scaled (2:1 ratio) 
by pixel replication into the new L 0 , the initialization for the next resolution k - 1 layer. 

14 The resolution layer is updated k = k - 1 , and the algorithm proceeds by going again to 

15 block 125. If k = 1, the result L T i, is the final output of the algorithm 120. 

1 6 Figure 7 is a block diagram of an alternative embodiment of an algorithm routine, 

17 denoted as 120a. The routine 120a may be applied to solve the equation: 

18 u" +l =u" +r\aL n -D") 



ubject to the constraint u n t] e 3. 



HP 10005732 



1 The routine 1 20a begins with block 151, where the kernels D 2 and g, and the 

2 counter k are initialized. 

3 In block 153, the routine 120a calculates: 

4 D" =g*g*(u n -w 0 ). 

5 In block 155, the routine 120a calculates: 

6 L" = D 2 *(u" -u 0 ). 

7 In block 1 57, the routine 120a performs the gradient steps to solve u" } +x as noted 

8 above. 

9 In block 158, the routine 120a determines if k>l, and if so, processing moves to 

10 block 159, where k is set to k-1. Processing then returns to block 153. Otherwise, the 

1 1 routine 120a ends. 

1 2 The penalty function as shown in Equation (5) tends to create halos in the 



13 resulting image. Figure 8 explains the origin of those halos through a one dimensional 

14 example. Figure 8 shows an original signal 160 that is outside of gamut 162 (delimited in 

15 the gray range by dotted lines 163 and 164). Projecting the signal 160 onto the gamut 162 

16 will result in a constant value, denoted as "High" (163) and loss of all detail. Figure 8 

17 also shows the result of processing the original signal 160. Dashed line 165 represents 

18 the result of scaling the signal 160 into the allowed range of the gamut 162. All the 

1 9 details are preserved, but with a smaller contrast. As opposed to point operations, the 

20 space dependent approach represented by Equation (5) yields a signal 1 66 (the solid line) 

21 that preserves the details with high contrast. However, halos occur near strong edges 167, 

22 which means that near the edges 167 there is a slow transition from low to high values. 

23 In order to avoid these phenomena, the penalty term may be modified, and robust 

24 estimation tools may be used. The original penalty term in Equation (5) may be replaced 

25 by: 

E(u)= {( A (D)+a/> 2 {VD|)Vfi. (8) 

n 

26 which for p l (x)= p 2 {x) = x 2 coincides with Equation (5). If the function p(x) grows 

27 slower than x 2 as x— > <x>, behavior near strong edges improves. Good candidates for p{x) 

28 are /?(x) = |x| or p(x)= ^jl + x 2 . 
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1 A different and simpler (linear) approach with similar robust behavior involves 

2 solving the original Equation (5) twice, with two different values of a. A solution with a 

3 small a may be denoted as u sma u and a solution that corresponds to the high value of a 

4 may be denoted as Uh lg h- The solution u sma ii has smaller contrast at areas with small 

5 details, yet has almost no halos. On the other hand, the solution Uh, g h preserves the small 

6 details, but at the expense of strong halo effects. By averaging these two results (u sma ii 

7 and Uhigh) in a spatially adaptive way provides a simple, improved solution. The 

8 improved solution is therefore: 

U final A = A k > j\ 'small [K /ft " >*fc / JKg* [*, /'] 

9 The weight w[k, j] should be close to one near strong edges, and close to zero in relatively 

10 smooth regions. In an embodiment, 

w [ kJ ] = 1 

\ + /3\Vg*u 0 \ 2 

1 1 provides a robust estimation. 

12 Halo problems have been recently dealt with in relation to Dynamic Range 

1 3 Compression. Solutions proposed included anisotropic diffusion and robust filtering. 

14 The halo-related solutions described herein are solutions to the same halo problem. 

15 Figure 9 is a flowchart showing an operation 180 of the gamut mapper 105. The 

16 process begins in start block 185. In converter block 195, the raw image S is converted 

1 7 from electrical signals having arbitrary values (colorimetric values) to locally defined 

1 8 color descriptors. 

19 In block 200, the converted image signal is decimated to form an image pyramid. 

20 A 2: 1 decimation scheme may be used, for example. Subsequent steps in the operation 

21 180 begin with the coarsest resolution layer of the image pyramid using a space varying 

22 algorithm. In block 205, the gamut mapper 105, using a space varying algorithm such as 

23 that represented by Equation (5) calculates the image deviation for the coarsest resolution 

24 layer. This process may involve calculating a first variation (block 210), determining a 

25 gradient descent flow (block 215), and solving the resulting gradient subject to a 

26 constraint (block 220). In block 225, a determination is made if the current resolution 

27 layer is the last (finest) resolution layer. If additional layers remain for processing, the 

28 process 1 80 returns to block 205. Otherwise, the process 1 80 moves to block 230, where 
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1 local color descriptors are converted to the required output format. In block 245 the 

2 process 180 ends. 

3 The above-described space-varying algorithm 120, as represented by Equation (5), 

4 for example may be executed using a general purpose computer system 250 as shown in 

5 Figure 10. The computer system 250 includes an input device 260, a computer unit 270, 

6 a display 280, a keyboard 290, and a storage device 300. The storage device 300 may be 

7 a CD-ROM drive, for example. Program code instructions for executing the algorithms 

8 120, 120a and 180 may be stored on a CD-ROM 305, or other computer readable 

9 memory, which is read by the storage device 300. The program code read out of the CD- 

10 ROM 305 are executed by a CPU of the computer unit 270. The CPU may perform the 

1 1 processing shown by the flowchart of Figure 9. The algorithms 120, 120a, and 180 also 

12 may be performed in a camera or printer hardware. 
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